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Abstract. We study collective behavior of Fermi gases trapped in various external 
potentials, including optical lattices (OLs), in the framework of the mean-field 
(hydrodynamic) description. Using the variational method, we derive effective 
dynamical equations for the one- and two-dimensional (ID and 2D) settings from the 
general 3D mean-field equation. The respective confinement is provided by trapping 
potentials with the cylindrical and planar symmetry, respectively. The resulting 
equations are nonpolynomial Schrodinger equations (NPSEs) coupled to equations for 
the local transverse size of the trapped states. Numerical simulations demonstrate 
close agreement of results produced by the underlying 3D equation and the effective 
low-dimensional ones. We consider the ground state in these settings. In particular, 
analytical solutions are obtained for the effectively 2D non-interacting Fermi gas. 
Differences between the ID and 2D configurations are highlighted. Finally, we analyze 
the dependence of the ID and 2D density patterns of the trapped gas, in the presence 
of the OL, on the strengths of the confining and OL potentials, and on the scattering 
length which determines the strength of interactions between non-identical fcrmions. 



PACS numbers: 03.75.Ss, 31.15.E, 74.20.Fg 
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1. Introduction 

In recent years, ultracold atomic gases have been explored as the means of experimental 
and theoretical simulations of various aspects of many-body physics [1]- [5]. In particular, 
conventional solid-state physics can be emulated by gases of fermionic atoms. The 
possibility to adjust the depth of the optical lattices (OLs), into which the gas is usually 
loaded, and the strength of interactions between atoms via the Feshbach resonance, 
opens the way to a detailed experimental study of these settings [6]- [7J. Such possibilities 
have inspired, in particular, many works dealing with the BCS-BEC crossover [8]-|llj. 
which occurs on the route from weak to strong attraction, being particularly relevant as 
a simulator of superconductivity. 

In these contexts, the density-functional theory (mean field) and the local-density 
approximation have been used with great success for the description of quantum gases 
close to their ground states [12]- [IB], resulting in effective time-dependent nonlinear 
Schrodinger equations. In the framework of this approach, hydrodynamic equations 
for the superfluid Fermi gas at zero temperature were derived in Ref. [19], and, in a 
similar way, the hydrodynamic equations which remain valid in the presence of temporal 
variations of control parameters were derived in Ref. [20J. 

The possibility of highly confined systems in one and two dimensions has stimulated 
the search for dimensionally reduced descriptions. In particular, quasi-one-dimensional 
Fermi systems described by many-body Hamiltonians have been analyzed, following 
seminal works [2T1 [22] and a more recent one [23]. A variational approach for ID and 
2D Fermi gases in the presence of OL potentials was developed in Ref. [21], and a similar 
technique for ID and 2D Fermi-Fermi mixtures was proposed still earlier in works [25] 
and [26] . 

In this paper we address Fermi gases in the BCS regime in nearly 2D and ID 
configurations, i.e., disk- and cigar-shaped ones, respectively. The configurations 
may also include inner OL potentials. To reduce the 3D nonpolynomial Schrodinger 
equation (NPSE) for the Fermi gas, known from the hydrodynamic description, to the 
corresponding 2D and ID equations, we use the variational method, similar to that 
employed in Refs. [27] 128] to reduce the 3D Gross-Pitaevskii equation for the Bose- 
Einstein condensate to the respective ID and 2D bosonic NPSEs. Unlike variational 
methods usually employed in the description of confined fermion systems, which 
postulate a constant width of the wave-function profile in the confined dimensions, the 
present method accounts for variations of the confinement width, resulting in significant 
corrections to the density profile. We demonstrate that our improved approximation is 
essentially closer to results of the 3D integration. 

As a starting point, we derive an effective hydrodynamics equation for a 
multicomponent Fermi gas with balanced populations, generalizing those found in the 
literature for a gas with two balanced populations, considering the BCS limit the density 
of energy. Next, we derive a effective equations for the disk-shaped trap, in the form of 
the respective equation for the fermionic mean-field wave function coupled to an equation 
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for the transverse width of the trapped state. Analytical solutions of the 2D equations 
can be found for the non-interacting gas; in the presence of the interactions between 
fermions with spins up and down, we develop numerical solutions, and demonstrate 
that the reduced equations reproduce the results obtained from the full 3D equation 
with a very good accuracy. We further extend the analysis to include the 2D (in-plane) 
lattice and superlattice potentials. Then, we examine the reduction of the 3D NPSE to 
ID equations in the case of the cigar-shaped confinement. The quasi-lD configuration 
including the axial OL potential is investigated too. 

2. The model 

2.1. Effective hydrodynamic equation for two balanced populations 

We consider a weakly interacting dilute Fermi gas of N atoms of mass m and atomic 
spin s, which form the BCS superfluid, loaded into an optical trap. The gas is taken 
at zero temperature, with balanced populations of the 2s + 1 possible spin states, as 
specified in detail below. To do this, we start with the case of balanced populations, for 
which an effective hydrodynamic equation was derived from the time- dependent density- 
functional theory for the Fermi gas, in the regime of the BCS-BEC crossover, by Kim 
and Zubarev in Ref. [T9] : 



is the chemical potential, related to the energy per particle e (n). In Ref. [19] it was also 
shown, for elongated cigar-shaped potential traps, that the von Weizsacker's correction 
to the Thomas-Fermi kinetic-energy density (which corresponds to the hydrodynamic 
approximation) [29], £tf (n) = (3/10) nh 2 (3n 2 n) 2 ^ 3 / m, may be important even for the 
gas with a large density of fermions, |Vn| /n 4 / 3 <C 1 (this is the applicability condition 
for the hydrodynamic approximation). Such an improved hydrodynamic approximation 
corresponds to the first two terms of the Kirzhnitz's gradient expansion [30]. For large 
smoothly varying densities in the free space, the respective correction is negligible; 
however, for strongly confined systems the additional term (the kinetic energy) may 
be comparable to the main one, hence it cannot be neglected, as shown below by our 
variational approach. 

A problem also arises at edges of the trapped superfluid, where |Vn| /n 4 / 3 ceases to 
be small. To overcome this problem, one may consider the trapped gas with large number 
of particles, lessening the relative significance of the edges. The case of a relatively small 
number of particles (N ~ 10 2 ) is considered below too, just with the aim to test the 
accuracy of the variational method. 




(1) 
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2.2. The BCS regime 

For a fermionic gas composed of two spin states with balanced populations and a 
negative scattering length of the s-wave collisions between the particles in different spin 
states, a s < 0, the BCS limit corresponds to fcp|a s | <C 1, with the Fermi wavenumber 
k F = {3n 2 n) 1 / 3 . In this regime, the energy per particle, e, is represented by the expansion 
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where e F = h 2 k F / (2m) is the Fermi energy. From Eqs. (j3J) and (|2j) it follows that the 



chemical potential takes the form 
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where the first term corresponds to the effective Pauli repulsion, and the following 
ones to the superfluidity due to collisions between the fermions in different spin states. 
Substituting the latter expression in Eq.([T]), and keeping only the first collisional term, 
we obtain the known nonlinear Schrodinger equation for the fermionic superfluid (see, 
e.g., Refs. [iHlOE]), 

2wh 2 a 
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(5) 



where the last term is similar to that in the Gross-Pitaevskii equation for bosons, but 
with an extra factor of 1/2, as the Pauli exclusion principle allows only atoms in different 
spin states interact via the scattering. Equation (jSJ) implies equal particle densities and 
phases of the wave functions associated with both spin states. 



2.3. The effective hydrodynamic equation for multiple balanced populations 

We now aim to extrapolate the previous case to systems with multiple atomic spin 
states, <7j, associated with 2s + 1 values of the vertical projection of the spin. To this 
end, we treat the atoms in each spin state (cTj) as a fully polarized Fermi gas, and include 
interactions between atoms in the different spin state (assuming a single value of the 
scattering length, a s ), which leads to the following system of equations: 
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where ^fj is the wave function of the superfluid associated with spin projection o~j, 
such that rij (r, t) = \^fj (r, t) | 2 is the respective particle density, and V(r) an external 
potential, which is assumed to be identical for all the spin states. A relevant example 
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of a higher nuclear spin, s = 9/2 (decoupled from the electron states), in a degenerate 
Fermi gas of 87 Sr, was demonstrated in Ref. [32] . 

In the case of fully locally balanced populations, the density of particles is the same 
in each component, n x = n 2 = ... = ri 2s +i, hence the total density is n = (2s + l)nj. 
Assuming also equal phases of the wave-function components, we define a single wave 
function, \1/ 



y/2s~+T 1 $ j , hence Eq. 
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(7) 



with density n 3 o (r, t) = |\& (r,t)| 2 and scattering coefficient 

g = 8sivh 2 a s /(2s + l)m. (8) 

In particular, the fully polarized gas, with the interactions between identical fermions 
suppressed by the Pauli principle, formally corresponds to s = (hence, g = 0, as given 
by Eq. (JS|)). Unless specified otherwise, we set s — 1/2 below, focusing on the case of 
balanced populations in two components. 



2.4- The Lagrangian density 

Equation can be derived, as the Euler-Lagrange equation 
5C dC d dC _ dC 




= (9) 
Sq* d$/* dtd(d t ^>*) "<9(W*) ' v; 

from the corresponding action, S = f dtdrC, with the Lagrangian density 
£ = i'n ( - ) - |^|V*| 2 - V(r)n 3D - 

2m5\2s + lJ - \ gnl ^ (10) 

where the asterisk stands for the complex conjugate. Similar Lagrangian formalisms have 
been used, in the context of the density-functional theory, in diverse settings [131 12S1 [33] . 

Our goal is to reduce the 3D description to one and two dimensions, using the 
variational approach. As said above, this derivation may be compared to that for the low- 
dimensional NPSE from the 3D Gross-Pitaevskii equation for the bosonic gas [27j [28] . 
The approach is based on the factorization of the 3D wave function, substituting the 
corresponding ansatz into the action, and integrating out the factor corresponding to 
the direction(s) orthogonal to the confined dimensions. 



3. The two-dimensional reduction 

3.1. The derivation of the 2D equations 

Here we aim to derive effective 2D equations for the Fermi gas in a disk-shaped trap. For 
this purpose, we consider the external potential composed of two terms: the parabolic 
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(harmonic-oscillator) one accounting for the confinement in the z direction, transverse 
to the disk's plane, and the in-plane potential, V 2 n: 



V(r) 



-mu) 2 z z 2 +V 2 D (x, y) 



As said above, the initial ansatz assumes the factorization of the 3D wave function into 
a product of functions of z and (x,y), the former one being the ground state of the 
harmonic-oscillator potential [T8] : 

V(r,t) = 1 exp | J]d>(x,y,t), (12) 

subject to the unitary normalization, with transverse width £(z,t) considered as a 
variational function, while the 2D wave function, <p, is normalized to the number of 
atoms. Therefore, the reduction from 3D to 2D implies that the system of equations 
should be derived for the pair of functions <j>(x,y,t) and £ (x,y,t), using the reduced 
action functional to be derived by integrating the 3D action over the ^-coordinate (cf. 
the derivation for the bosonic gas [271 EH]): 

^2D = J dtdxdyC 2 v, (13) 
where the respective Lagrangian density is 
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with 2D density n 2D (x,y) = |0(x,y)| 2 , C 2D = (3/5) 1 / 2 (6/(2s + l)) 2 / 3 vr. The last two 
terms in Eq. ( 1T4"|) result from the reduction to 2D: the first among them shows that 
stronger confinement implies a higher energetic cost, while the quadratic dependence on 
£ in the last term originates from the original confining potential acting in the z direction. 
These two terms are relevant for configurations with an inhomogeneous spatial density, 
while for spatially homogeneous solutions they may be omitted. 



Atoms 


u c [Rz] 

a c =0.1xl0- 6 m 


w c [Hz] 

a c =1.0xl0 -8 m 


w c [Hz] 

a c =10.0xl0- 6 m 


6 Li 


10.5 x 10 5 


10.5 x 10 3 


105. 


40 K 


1.58 x 10 5 


1.58 x 10 3 


15.8 


84 Sr 


0.75 x 10 5 


0.75 x 10 3 


7.5 



Table 1. Relevant values of the characteristic frequency, u> c , for 6 Li, 40 K and 84 Sr, 
assuming three possible values of the characteristic length, a c . 



To present the effective 2D equations in a dimensionless form, we define a 
characteristic length, a c , and a characteristic frequency, u c , related through the atomic 
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mass (m), u c = Ti/{ma 2 c ). In this way, the physical units of the system depend on 
fix the value of a c (or u c ) for particular atomic species, see Table [TJ Thus, we rescale 
the variables and constants as t = u c t/2, {X, Y} = {x, y} /a c , d = a c 0, p 2 n = |^| 2 , 
C = £/a c , U 2 t> = 2V r 2D /(^w c ), G 2 b = 2 7/2 v/7r 'sa s / [(2s + l)a c ], and w z = u z /u c . In this 
notation, time t is expressed in units of 2/u c , and 2D coordinates {x, y} in units of a c . 
Then, Eqs. f JT3|) and ffbil) are transformed into 



S = h / drdXdY£ 2D (15) 



and 

£ 2D = \ (^d T d - dd T d*) - | W| 2 - U 2D p 2D 

^T^CWsb + -7^P2D + ^P2D + -W,C P2D , (16) 

where the part in the square brackets corresponds to the local energy density, e 2 d, when 
the system is in a spatially homogeneous state in the absence of external potentials: 

3 ri 5/3 . C 2 D 2 1 1 2/2 /i>7\ 

^2D = ^^C , 2DP 2 D + T^P2D + ^2"P2D + ~^ Z C p2B ■ (17) 

The latter equation is a relation between the energy and particle densities, demonstrating 
how different interaction terms affect the local stability of the matter wave. 

The Euler-Lagrange equations produced by varying action S with respect to d and 
C, take the form of 



-Vi + ^2 D + + + ^ + ^ 



0, 



wlC* - \c 2T > |tf| 4/3 C 4/3 - l -G m -1 = 0, (19) 



cf. Eq. (J9J). Algebraic equation ffl9|) for ( cannot be solved analytically, therefore we 
used the Newton method to solve it numerically The necessity to find ( at each step 
of the integration is a numerical complication of a minimal cost compared to the 3D 
integration of the underlying equation ([7j). 

Nevertheless, for the non- interacting gas, with G2D = 0, Eq. ( Fl9l) is a cubic equation 
for C 4 ^ 3 , which admits an analytical solution, the single physically relevant real root being 

1/3 + 2 1 / 3 



-, 2/3 C 4/3 (fi) = —^=^-T7, + ^ %m ! , (20) 



where fi = 2C 2D |^| 4/3 /(5ro 2 / 3 ). This root is real at p 2D /w z < 5 9 / 4 (2s + 
l)/(2 7//2 3 1//4 7r 3 / 2 ) Ri 0.45(2s + 1), i.e., for sufficiently low densities and the strong 
transverse confinement. Further, in the low-density limit, the approximation for the 
root is 

1 / Con , . , a in \ , _ . . 
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Figure 1. (Color online) The energy density and transverse size, £2D and £, versus 
the 2D atomic density, p 2 D- (a) and (b): The dependences are shown for five different 
values of the spin, s = 0, 1/2, 3/2, 5/2 and 7/2, at the fixed negative scattering length, 
Q>s/°>c = —0.2, and fixed confinement strength, vo z = 5. (c) and (d): The same 
for five different values of the scattering length, a s /a c = 0, —0.05, —0.10, —0.15 and 
—0.20, while s = 1/2 and w z = 5 are fixed, (c) and (f): For five different values 
of the transverse-confinement strength, w z = 1,2,3,4 and 5, at fixed s = 1/2 and 
a s /a c = —0.2. The optical lattice is absent in all the cases. 



the substitution of which into Eq. ( fl8|) results in the following NPSE for d: 



id T $ 



4/3 ^2D| 



15^ 2 2/3 



tf, (22) 



where the chemical potential is shifted by a constant term, 2w z . In the limit of —> 0, 
Eq. (|2ip yields C = -o?^ 1 / 2 , which corresponds to the width of the ground state of the 
harmonic oscillator. 

The plots in Fig. [1] display the dependence of the energy density e 2 D and transverse 
width C on the 2D density p 2 D for the states uniform in the (x, y) plane, as produced 
by Eqs. (IT7|) and Eq. ffl9l) . respectively (recall that p 2 D = |^| 2 )- In particular, Figs. [T](a) 
and H](b) show how these dependences are affected by the fermion's spin, s. In Fig. Ufa), 
the energy density of the polarized gas, with s = 0, is a growing function of the atomic 
density, since in this case Eq. (JEJ) gives g — 0, i.e., the Pauli exclusion principle suppresses 
the interaction between the fermions, as mentioned above. When the gas is composed 
of equal numbers of atoms with spins up and down (s = 1/2), the attraction (a s < 0) 
between atoms with opposite directions of the spin causes a decrease in the energy 
density, in comparison to the fully polarized gas. For the cases of s = 3/2, 5/2 and 7/2, 
it is observed that a maximum of the energy density is attained at some atomic density 
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Figure 2. (Color online) (a) The maximum of the 2D energy density, £2D,max! a s a 
function of the negative scattering length, a s /a c . (b) The 2D atomic density, p2D,max, 
at which £2D attains its maximum. The curves in both plots correspond to five different 
values of spin s, as indicated in the inset in (b). The confinement strength is fixed to 
be va z = 5, and the optical lattice is absent here. 



(p2D,max), due to the greater contribution of the attractive interactions in comparison 
with the quantum pressure. In agreement with this trend, p2D,max is seen to be smaller 
for higher values of s, as the greater number of possible spin states allows each atom 
to directly interact with a larger number of atoms in other spin states. Furthermore, 
Fig. QJb) shows that (for the same parameters as in Fig. d^a)) the gas with s — 1/2 
slowly shrinks in the transverse direction (( gradually decreases) with the increase of the 
2D density, after attaining a weakly pronounced maximum at certain value of /?2D- It is 
also observed that, at the vanishingly small density, the transverse size does not depend 
of spin s, in agreement with Eq. ( )2lT) . On the other hand, Figs. [I^c) and [Hd) show 
that, for s = 1/2, the energy density remains a growing function of the atomic density 
for all the values of a s /a c considered, and, naturally, the growth is steeper for smaller 
absolute values of the negative scattering length. Further, the curves in Fig. []Je) show 
that the energy density grows steeper when the confinement strength, w z , is higher, and, 
naturally, the tighter confinement leads to the smaller transverse size, as seen in Fig.Q^f). 
In summary, for s > 1/2 and in the range of the particle density (p2d) considered here 
(see Fig. [l^a)), the energy density (£2d) is an increasing function of the particle density 
up to a certain value of the density (p = P2D,max) ; where e^d attains a maximum, which 
it followed by the decrease of the energy density. It is also observed that attaining 
this maximum is always preceded by the gradual tightening of the gas in the confined 
direction and increase of the particle density. In dynamical simulations, the system 
suffers the collapse if the initial particle density exceeds the maximum value p2D,max 
(in a nonuniform configurations, the particles are pulled towards the point where this 
maximum is reached), i.e., p2D,max is a separatrix bordering between stable and unstable 
uniform states, cf. similar results reported in Refs. [261 GUI US] . 

Figure [2] shows the maximum value of the energy density as a function of the 
strength of the attraction between the fermions in the different spin states. It is 
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Figure 3. (Color online) The 2D radial density, P2T>(r), as obtained from the full 3D 
equation, and from the 2D reduction based on the variational approximation ("VA") 
for the weak isotropic harmonic-oscillator trap acting in the 2D plane (w x = vo y =0.1, 
zu z = 2), and N = 1000. Different curves correspond to the indicated values of 
a s /a c = 0,-0.05 and —0.1. The optical-lattice potential is absent in this case. 



seen that the enhancement of the attraction causes a transition to lower energy and 
atomic densities. This behavior is stronger pronounced for lower values of the spin, 
as a result of the larger relative strength of the Pauli repulsion, leading to greater 
variations in the energy- density maximum. It should be stressed, though, that these 
results are approximate due to limitations imposed by the BCS regime (overall, the 
weak interaction). Actually, most calculations reported in this work were performed for 
the energy densities far from the maximum value. 

After the consideration of the spatially uniform states, the main objective of our 
analysis is to identify the ground state of the fermion gas loaded into the 2D potential 
consisting of the harmonic-oscillator trap and square OL, 



where K x y = 2tt/\ x ^ are the lattice wavenumbers. We start the presentation of the 
results for nearly-2D trapped Fermi gas in the absence of the OL. 

3.2. The trapped quasi- 2D gas in the absence of the optical lattice 

First, we have verified the accuracy of the 2D reduction by comparing results generated 
by this approximation versus those obtained by integrating the 3D equation (0). The 
ground state was found by means of the imaginary-time integration based on the fourth- 
order Runge-Kutta algorithm. The comparison is produced in Fig. [HI where the radial- 
density profiles are plotted for the isotropic 2D harmonic-oscillator trap [w x = w y ), 
showing an excellent agreement between the 2D and full 3D descriptions. Actually, this 
is one of main results of the present work. 



U 2D = ™ 2 X X 2 + w 2 Y 2 + A x sm 2 (k x X) + A y sin 2 (n y Y) , 



(23) 
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Figure 4. (Color online) Results for the quasi-2D gas trapped in the combination of 
the optical lattice and weak harmonic-oscillator potential. The plots display density 
profiles, p2u(x 1 y), for four cases: (a) a s /a c = and w z — 2, (b) a s /a c = —0.2 and 
137, = 2, (c) a s /a c = and w z = 8, and (d) a s /a c = —0.2 and w z = 8. The remaining 



parameters are N = 1000, vj x = w y = 0.1, A x = A y = 10, and = 



9. 



3. 3. Effects of the optical lattice 

The plots in Fig. H] show the density distribution in the ground state of the gas loaded 
into the 2D harmonic-oscillator potential combined with the 2D square OL (see Eq. 123]) . 
Naturally, density peaks coincide with local minima of the OL potential. Figures H](a) 
and H](b) demonstrate that the attractive interaction leads to a decrease in the number 
of peaks visible in the density profile, and to an increase in their heights. On the 
other hand, Figs. H](c) and 11(d) show that, for the same interaction strength as above, 
tighter transverse confinement results in a greater number of peaks, with a lower atomic 
density at them. Due to the symmetry of the OL, the number of visible peaks must be 
in multiples of 4. We define as "observable" an area with density p 2 D > 0.1. The plot 
of Fig. [5(a) shows the dependence of the observable number of peaks on the transverse- 
confinement strength w z , for a s /a z = (green circles) and a s /a z = —0.2 (blue squares). 
At all values of w z , the number of visible peaks is higher in the a s = case, with 
a clear nonlinear dependence that is characterized by intervals in the scale of in w z 
corresponding to particular numbers of the peaks. This feature is particularly noticeable 
in the range of w z = [5, 10]. 
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Figure 5. (Color online) (a) The number of density peaks, located at minima of the 
OL potential, as a function of confinement strength w z . The area occupied by the gas 
is defined as that with p 2 D > 0.1. The blue and green lines correspond to a s /a c = —0.2 
and a s /a z = 0, respectively. Panels (b) and (c) display, severally, detailed pictures for 
a s /a c = —0.2 and a s /a c = 0, in seven different ranges of density/^D, as indicated 
in the inset. The parameters are N = 1000, w x = vj v = 0.1, A x = A y = 10, and 
A^ = A,, = 9. 
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Figure 6. (Color online) The same as in Fig. [4] for the supcrlattice built of two 
triangular OLs with angle 8 = 8° between them: (a) a s /a z = —0.2, w z = 2; (b) 
a-s/a-z = 0, w z = 5. The other parameters are N = 1000, vo x = w y — 0.1, and At = 5, 
A T = 10. 



Figures \Mb) and Efc) show detail of the peaks characterization, which amounted 
to counting the number of peaks in seven specific zones of values of />2D ; [0. 1,1 [, [1,2 [, 
[2,3[, [3,4[, [4,5[, [5,6 [ and [6,7[. The findings uncover intricate changes in the density 
patterns, following variations of the system's parameters. In particular, Fig. [5]^b) shows 
that, for a s /a z = —0.2, distinct zones of the values of w z may cluster at the same 
number of peaks, suggesting that a given number of peaks may be realized in various 
configurations. The plot in Fig. [5^c) shows the non-interacting case, a s = 0, with no 
visible peak in the last two ranges, which implies greater expansion of the gas in the 
2D plane. In each range, the variation of parameters significantly alters the distribution 
of peaks. The multitude of the different coexisting robust multi-peak patterns suggests 
that this setting has a potential for the use as a data-storage system. 



One- and two-dimensional reductions of the mean-field description of degenerate Fermi gasesl3 



Figure [6] displays the results for the Fermi gas in the presence of a superlattice 
formed by a superposition of two triangular lattice, each formed by the set of three 
pairs of counter-propagating laser-beams, with angles 60° between them, and equal 
amplitudes and wavelengths. To build the superlattice, we fix the first pair of beams 
in the x-direction in both triangular lattices, and then rotate the two lattices clockwise 
and counterclockwise by ±4°. The overall shape of the thus formed superlattice remains 
hexagonal. It is observed that the density maxima are most visible within the central 
hexagon. The set of different patterns trapped in the superlattice is richer than in the 
case of the square OL considered above, as different potential minima of the superlattice 
potential have different depths. 



4. The one-dimensional reduction 

4-1. The derivation of the ID equations 

To perform the reduction of the mean-field description from 3D to ID, we consider an 
external potential formed by the harmonic-oscillator potential in the (x, y)-plane, which 
represents the cylindrically symmetric magnetic confinement, and an axial (z- dependent) 
potential: 

V(r,t) = ~mu 2 t r 2 + V m (z,t), (24) 



where r = \fx 2 + y 2 . Again, the initial ansatz assumes the factorization of the 3D wave 
function into a product of r- and ^-dependent functions. The former one is taken as the 
ground state of the 2D harmonic oscillator, hence the ansatz is 

*< r - i )=7^(M) exp ("^b) /(z ' i) ' (25) 

where the ID wave function, / (z, t), is normalized to the number of atoms, N, and the 
axial density is defined as nm = \ f\ 2 - Variational functions / and a can be determined 
by the minimization of the 3D action after performing the integration in the (x, y) plane, 
cf. the derivation of the ID NPSE for the bosonic gas [27]. This procedure gives rise to 
the effective ID Lagrangian density, 

Cm = i\ (fdtf - fd t P) - ^\d z f\ 2 - V m n m 

3 5/3 g n ^ 2 1 2 2 

_^^ 1DniD + W" 1D + 2W" 1D + 2^ niD \ ' (26) 

where Cm = (3/5)(67r/(2s + 1)) 2 ^ 3 , s is the fermionic spin, as before, and the last two 
terms are the result of integration, which keeps the 3D characteristics of the underlying 
system: When the gas is homogeneous and the external potential is absent, the latter 
term is irrelevant. 

Similar to the 2D case, the normalization is performed with characteristic length 
scale a c and frequency co c , so that Z = z/a c , ip = a x J 2 f ■, U\v> = 2Vid/(^c)> 
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Gid = 8sa s /(2s + l)a c , 77 = o~/a c , and w t = oo t /uj c . In addition, we define the 
dimensionless ID atom density, Pid = |^| 2 - Then, the normalized Lagrangian density is 



C 



ID 



^ ~ 5/3 , 1 ^2 , „2^,2 



CidPid + tt^^idPid + w t ?7 Pid + — PlD 



5r/ 4/3 ""-^ ' 2 p 

The ID energy density £id for uniform states is 

3^5/3.1^,2 22 1 

£lD ~ 5rf/3 1DplD + 2^f GlDplD + OTt?7 PlD + ^2 plD > 



(27) 



(2? 



cf. Eq. (I17p in the 2D case, which makes it possible to compare contributions of the 
different interactions involved. In fact, the comparison is not straightforward, because of 
the dependence of rj on pxo\ however, for very low densities 77 may be assumed constant, 
in which case the last two terms dominate in em. Now, varying the action with respect 
to ip, we derive the respective Euler-Lagrange equation (similar to the Eq. [9]): 



id T ip 



c 



ID 



ID 



^4/3 nqi. 

which is the ID equation of motion, with powers of the transverse width, 77, in the 
collisional and Pauli terms being different from their counterparts in the 2D equation, cf. 
Eq. (j!8p . Further, varying the action with respect to 77 gives rise to the corresponding 
Euler-Lagrange equation, 



1 4/3 



G 



ID 



1 2 2 

— + &tV 

7] z 



(29) 



2 4 
w t rj 



14/3^2/3 



0. 



(30) 



which establishes the dependence of the transverse size, 77, on density piD- Generally, 
this equation should be solved numerically. In the low-density limit, it has a trivial 
solution, if = l/w t , which, if substituted in Eq. (1291) allows one to derive a closed- form 
equation for ip [similar to the situation for the same limit in the 2D case, cf. Eq. (122]) ]: 



id T ip 



-dl 



U 1D + C 1D w 2 t /3 \ij\ 4/3 + w t G 1D 



2zu. 



1>, 



(31) 



with the lowest-order nonlinear term accounting for the Pauli repulsion of the same type 
as in the 2D equation ( J22l ). 

The plots in Fig. [7] show the dependence of the energy density , Em, and transverse 
width, 77, on the ID density, pin, as obtained from Eqs. ( |28l) and ( |30l) . respectively. The 
dependences are similar to their 2D counterparts shown in Fig.HJ i.e., the energy density 
increases with the atomic density, the increase being stronger at higher values of ta t , 
and with the strength of the repulsive interactions (a s > 0). The similarity is also true 
for the dependence of the transverse width, 77, if compared to that for width, (, in the 
2D model. However, estimating the 3D density, 77,30, in the case of piD ~ P2D (which 
also implies that ( ~ 77), we find that the 3D density, as produced by the ID model, is 
larger by a factor ~ l/a c . Therefore, the range of values of 7730 corresponding to Fig. [7] 
is much broader than in Fig. [T], implying that, using the strong 2D confinement, the 
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Figure 7. (a) em vs. pm and (b) 77 VS. pir>. In both plots the five curves correspond 
to a s /a c = —0.2, —0.1, 0, 0.1 and 0.2, the confinement parameter being tut = 5. (c) eid 
vs. piD and (d) 77 vs. piD- In both plots the five curves correspond to w% = 1,2,3,4 
and 5, and the scattering length is fixed as a s /a c = —0.1. 



unitary regime can be achieved at much lower densities than with the ID confinement, 
which is explained by the fact that the Pauli-repulsion term is more important in the 
latter case. 

Again, we assume that the gas is subject to the action of the combination of the 
harmonic trap and OL in the axial direction, 

U m = w 2 z {Z - Z f + A z sm 2 {k z Z) , (32) 

where k z = In j\ z is the OL wavenumber, and A 2 /2 the OL period. The results obtained 
for such axial potentials are reported below. 

4-2. The trapped quasi-lD gas in the absence of the optical lattice 

As in the previous section, it is first necessary to check the accuracy of the 
approximation. To this end, in Fig. IH(a) we compare the ID density, pm, as produced 
by three different methods: the full 3D equation (|7j), the full ID approximation, and 
its low-density limit, see Eq. (|3T|) . Similar to the 2D case, we observe an excellent 
agreement between the 3D equation and the full ID approximation, both significantly 
disagreeing with the low-density limit. In particular, Fig. E^b) demonstrates that the 
approximation remains accurate in the presence of the interactions too. On the other 
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z z 



Figure 8. (Color online) The ID density profile, pm, as calculated through the full 
3D solution ("3D"), and by dint of the full ID reduction based on the variational 
approximate ("VA") (panel (a) also includes the limit case of the latter, corresponding 
to the"very- low-density variational approximation" ("VLDVA")). (a) Parameters are 
N = 100, a s /a c = 0, rot = 1 and w z = 0.1. (b) For three values of the scattering 
length, a s /a c = 0, —0.05 and —0.1, where the other parameters are N = 100, Wt = 2 
and w z = 0.1. (c) N = 500, and (d) N = 2000, for the other parameters a s /a c = —0.05, 
zut = 2 and w z = 0.1. 



hand, Figs. Ws) and (d) show that an increase in the number of atoms does not affect 
the agreement between the ID and 3D simulations. 

Figure[9^a) shows profiles of the ID density, Pid, for five different values of spin s, in 
the presence of the confining harmonic-oscillator potential acting in the axial direction 
(z). The smallest width of the profile is observed at largest s, which is a result of the 
relatively weak Pauli repulsion versus the attractive interactions, as a consequence of 
the lower spin degeneracy. Figure Mjo) shows axial expansion of the trapped Fermi gas 
with the increase of the strength of the transverse confinement (w t ). This trend is a 
consequence of the rapid enhancement of the Pauli repulsion. 

4-3. Effects of the axial optical lattice 

Typical ID-density patterns for the gas trapped in the OL are displayed in Fig. HUT a). 
Note that the increase of the confinement strength leads to attenuation of the central 
peak and generation of new side peaks, which is consistent with the elongation of the gas 
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Figure 9. (Color online) (a) The ID density profile. Five ID density profiles 
correspond to s = 0,1/2,3/2,5/2, 7/2 (as indicated in the inset) and Wt = 5. (b) 
The same, with the five curves corresponding to Tu t = 1,2,3,4, 5 (as per the inset) 
and s = 0. The parameters are TV = 100, a s /a c = and w z = 0.1. 
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Figure 10. (Color online) One-dimensional density profiles of the Fermi gas trapped 
in the OL with amplitude (depth) A z = 10 and double period X z = 20. (a) The five 
curves correspond to the confinement strength tut = 1,3,5,7,9 (as indicated in the 
top right inset). The parameters are N = 100, G = 0, and vo z = 0.1. (b) The five 
curves correspond to scattering lengths a s /a c = 0.2,0.1,0,-0.1,-02, as indicated in 
the panel. The parameters are N = 100, tu t = 5, and w z = 0.1. 

observed in Fig. Mjo). The lowering of the central peak is accompanied by the increase 
of the density between peaks, resulting in the reduction of their visibility. A similar 
behavior is observed in Fig. fTOT b). where the repulsive inter-atomic interaction makes the 
trapped state elongated along z, resembling the elongation induced by the strengthening 
transverse confinement, and vice versa for the attractive interaction. Figure [11] shows 
the strong dependence of the ID density profile piD on the OL's (double) period X z , 
for a fixed strength of the confinement potential. It is clear that the number of peaks 
in the ID pattern is larger for smaller X z . The increase of X z entails disappearance of 
peripheral peaks at particular values of X z . The present ID description makes it possible 
to identify these critical values. 

To further analyze the ground state of the Fermi gas in presence of the OL, we 
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Figure 11. (Color online) Onc-dimcnsional profiles of the Fermi gas trapped in the 
axial optical lattice with amplitude A z = 10. The plot shows density /?id versus the 
double lattice's period X z . The parameters are N = 100, G = 0, zut = 5, vo z = 0.1. 
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Figure 12. (Color online) Diagrams in the plane of the trapping frequencies, (tn t , w z ): 
(a) width A 2 of the area wherein /?id > 0.05; (b) the number of density peaks at the 
periphery of the trapped state, which are disconnected from the central peak. The 
parameters arc N = 100, a s /a z = —0.05, A z = 10, X z = 10. 



define the width of the trapped state, A 2 , as the size of the area with density Pid > 0.05. 
Figure [12] presents diagrams for characteristics of the ground state in the plane of the 
transverse and axial trapping frequencies, {m t ,m z ). The jumps observed in Fig. [T27 a) 
correspond to the emergence or disappearance of pairs of density peaks. Since the 
wavelength is X z = 10, the emergence or disappearance of the peak pair gives rise, 
respectively, to the increase or decrease of the width by A z = 10. Thus, the observed 
minimum width (A 2 = 20) corresponds to the profile with five peaks, and the maximum 
width (A z = 110) is represented by the profile with twenty three peaks. Note that, as 
seen in Fig. HOj the density between peaks may be very small. 

Further, we define that two peaks as being "connected" if the minimum density 
between them is greater than 0.05. Accordingly, Fig. [T2T b) shows the number of peaks 
"disconnected" from the central core. This characteristic is found by subtracting the 
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number of the connecting links between the peaks, defined as said above, from the total 
number of the observed peaks (excluding the central one). The diagram features four 
regions covered by similar colors, which correspond to four possible values of the number 
of isolated peaks, 0,2,4 and 6. These regions alternate with the variation of both w t 
and w z , demonstrating that the number of the isolated peaks may be maintained while 
varying the total number of the observed peaks. 

5. Conclusions 

We have demonstrated that the ID and 2D reductions of the 3D mean-field description 
of the degenerate Fermi gas, based on the Gaussian variational ansatz, provide very 
accurate results for the ground states of the gas in the disk- and cigar-shaped traps. The 
so derived reduced equations are useful, in particular, for modeling the BCS regime at 
low atomic density and large spatiotemporal variations of the external potential. For the 
ID case, the reduced equations provide results requiring the use of modest computational 
resources and allow one to quickly explore a vast volume of the parameter space. In the 
2D case, the necessary simulation time is still significantly lower than what is necessary 
for the 3D simulations. For the spatially uniform ground states in both dimensions, we 
have studied the width of the Fermi gas in the confinement dimensions, and the energy 
density as a function of the atomic density. In the 2D case, we have found that the 
energy density attains a maximum at a particular value of the atomic density, and the 
maximum shifts to lower values as the attractive interactions get stronger. While the 
analysis is approximate due to limitations of the model, it shows the correct dependence 
on the scattering length, confinement strength, and the phase-space degeneracy (i.e., 
the fermionic spin, s). 

For the gas trapped in the 2D OL, the analysis predicts the number of visible 
peaks in the density profile as a function of the scattering length and confinement 
strength. It was concluded that small variations of the latter parameter produce 
significant changes in the size and distribution of the peaks. In the ID setting including 
the axial OL, we have produced the phase diagrams for the Fermi gas, varying the 
confinement parameters. Discrete leaps in the diagrams correspond to the appearance 
and disappearance of a pair of peripheral peaks in the density profile. In fact, the 
number of peaks does not change in a wide range of the confinement parameters. 

In general, our model is able to include any type of the confinement and potential 
lattices in the ID and 2D settings, making the description of the ground state quite 
simple. In particular, we considered the example of the hexagonal superlattice built 
as a superposition of two triangular lattices. The latter setting may be essential for 
modeling configurations relevant to solid-state physics. The same approach may be 
further applied to superlattices built of hexagonal OLs, with the objective to emulate 
graphene-like superlattice. Beyond the study of the ground states, it may be interesting 
to build quasi-2D modes with embedded vorticity, cf . Ref . [21] . On the other hand, the 
reduced equations can be used to study the dynamics of the Fermi gas under the action 
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of external potentials slowly varying in time, similar to how it was done in BEC models. 
Finally, the same type of the description, including, if necessary, equations for the Bose 
gas, may be used to predict the ground state and analyze the low-energy dynamics of 
Fermi-Fermi and Bose-Fermi mixtures. 
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